!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! calculate x4, y4 from given nozzle wall function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! wall function is y=a+b*x+c*x^2+d*x^3
! slope of wall is dy/dx=tan(theta)=b+2*c*x+3*d*x^2
! slope of the slope of the wall is d^2 y/dx^2==2*c+6*d*x
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!   aw,bw,cw,dw    wall equation constant
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning:
!   all angle in rad, not degree
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine Wall34(lp)
   use VariableDef
   implicit none
   real*8::lp
   
   if (abs(thetae-thetaa) <= 0.01) then   ! line
      cw=0.0
      dw=0.0
      x4=(aw-y2+lp*x2)/(lp-bw)
      y4=aw+bw*x4+cw*x4**2+dw*x4**3
      theta4=atan(bw+2*cw*x4+6*dw*x4**2)
   else   ! curve
      x4=2*cw**3-9*bw*cw*dw+27*aw*dw**2+9*cw*dw*lp+27*dw**2*lp*x2-27*dw**2*y2
      x4=(x4+sqrt(4*(-cw**2+3*bw*dw-3*dw*lp)**3+x4**2))**(1/3)
      x4=-cw/(3*dw)+2**(1/3)*(-cw**2+3*bw*dw-3*dw*lp)/(3*dw*x4)+x4/(3*2**(1/3)*dw)
      y4=aw+bw*x4+cw*x4**2+dw*x4**3
      theta4=atan(bw+2*cw*x4+6*dw*x4**2)   ! slope angle
   end if
end subroutine Wall34
